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Abstract 

The cumulant representation of the Fourier path integral method is examined to determine the 
asymptotic convergence characteristics of the imaginary-time density matrix with respect to the 
number of path variables N included. It is proved that when the cumulant expansion is truncated at 
order p, the asymptotic convergence rate of the density matrix behaves like N~ ( - 2p+1 \ The complex 
algebra associated with the proof is simplified by introducing a diagrammatic representation of the 
contributing terms along with an associated linked-cluster theorem. The cumulant terms at each 
order are expanded in a series such that the the asymptotic convergence rate is maintained without 
the need to calculate the full cumulant at order p. Using this truncated expansion of each cumulant 
at order p, the numerical cost in developing Fourier path integral expressions having convergence 
order _/V~( 2p+1 ) is shown to be approximately linear in the number of required potential energy 
evaluations making the method promising for actual numerical implementation. 



I. INTRODUCTION 



The path integral methocP^ has proved to be an important and useful computational 
vehicle for obtaining thermodynamic properties of interacting many-particle systems in the 
quantum domain. In all path integral approaches the usual classical degrees of freedom in 
a system are augmented by an infinite set of path variables that effectively describe the 
quantum fluctuations about the classical trajectories. In actual simulations the infinite set 
of path variables is truncated, and an important issue is the rate of convergence to the exact 
quantum result as a function of the number of path variables actually included. 

In one approac h 12 * 4 * 7 ^ to path integration the imaginary-time propagator is discretized in 
coordinate representation using a large set of intermediate coordinate states along with 
the Trotter approximation.® The Trotter decomposition becomes increasingly accurate as 
the number of discretized path variables is increased. The asymptotic convergence rate of 
this discretized version of path integrals is knowrP to be 1/N 2 where N is the number of 
discretized points included. For problems where pair potentials are adequate, the conver- 
gence of the discretized method can be enhanced by, for example, using more accurate pair 
propagatorsP 

In this work we focus on an alternate path integral method 3 - 121 ^ where the quantum 
paths are expanded in a Fourier series, and the integration over all paths is replaced by a 
Riemann integral with respect to the Fourier coefficients. While exact results are obtained 
if the complete Fourier series containing an infinite set of terms is included, in practical 
applications, the number of coefficients included is truncated at N terms. In its primitive 
form this Fourier path integral method is knowrP^to converge asymptotically to exact results 
as 1/N. 

To enhance the asymptotic convergence rate of the Fourier method, a set of useful ap- 
proaches has been developed to include approximately the contributions from the coefficients 
excluded when the full Fourier series is truncated. The first of these methods has been named 
"partial averaging," 1 3 1 12 1 13 1 because integrals over the high-order Fourier path variables, the 
"tail integrals," are included in an average sense. Because partial averaging requires the 
evaluation of the Gaussian transform of the interaction potentials associated with a par- 
ticular problem, and because many interaction potentials used commonly in simulations 
do not have finite or readily available Gaussian transforms, alternative methods have also 
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been introduced that circumvent the need for a Gaussian transform. Among such methods 
are the gradient partial average methocP^ and the reweighted path integral techniqueP^ffiH 
The asymptotic convergence rates of the full partial average method^ and the reweighted 
method^ are both 1/iV 3 , whereas the gradient partial average method converges as 

In a previous publication^ we considered a suggestion by Singer 21 to fit the Lennard- 
Jones potential, which does not have a finite Gaussian transform, to a sum of two Gaussians, 
which have Gaussian transforms that are both analytic and finite. In that previous work we 
showed the fit potential to be an accurate representation of a one-dimensional Lennard- Jones 
system, and we explored the asymptotic convergence characteristics using a variety of Fourier 
path integral methods. In that work we found numerically, for the case studied, that the 
full partial average method reached its asymptotic limit more rapidly than the reweighted 
method. Recalling that the reweighted method has the same asymptotic convergence rate 
as the full partial average method, the results indicated that the full partial average method 
can be advantageous. 

The partial average method can be understood 313 ' as the first term in a cumulant 
expansio n 1 22 1 23 1 with respect to the tail series of the quantum density matrix. Motivated 
by the successful implementation of the partial average method using Gaussian fit poten- 
tials, in this work we explore the asymptotic convergence characteristics of the approach 
when cumulants are included to arbitrary rather than just first order. While the cumulant 
expansion for Fourier path integrals has been examined through second order in previous 
workED the properties of the full expansion have not been discussed before. The cumulant 
expansion by itself is not guaranteed to be convergent. In the current work we ignore poten- 
tial convergence issues with respect to cumulant order so that we can explore the asymptotic 
convergence characteristics of the cumulant expansion with respect to the number of included 
path variables N when the cumulant expansion itself is truncated at an arbitrary but finite 
order p. The analysis is complex, because at each cumulant order beyond p — 1, there are 
many terms having differing convergence characteristics with respect to N. We find many 
of the resulting terms cancel in a way that allows us to show that the resulting asymptotic 
convergence rate is 7V~( 2p+1 ). The analysis required is simplified by using a combination 
of algebraic and diagrammatic methods. The diagrammatic approach allows us to prove a 
linked-cluster theorem that tells us which of the large number of terms at each order in p 
survive cancellation. From the proved asymptotic convergence rate, the first-order cumulant 
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(partial averaging) converges asymptotically as l/N 3 , including the second-order cumulant 
enhances the asymptotic convergence rate to 1/iV 5 and so on. The derivation of the jV~( 2p+1 ) 
asymptotic convergence rate for cumulants of order p is a principal finding of the current 
work. 

The derivations of the asymptotic convergence rates rely on certain expansions at each 
cumulant order. We show that these expansions at each cumulant order, when truncated, 
have the same asymptotic convergence characteristics as the full cumulants. Furthermore, 
the truncated expansions at each cumulant order enable the development of path integral 
approaches that scale approximately linearly in the number of required potential energy eval- 
uations while retaining the N~^ 2p+1 ^ asymptotic convergence rate. This important finding 
implies the developments described in this work are potentially important from a numerical 
standpoint. 

The contents of the remainder of this paper are as follows. In the next section we present 
our theoretical developments. Included is a review of the Fourier path integral method, the 
full partial average method and the cumulant expansion. As examples we derive explicitly 
the asymptotic convergence rates of the first through third order cumulant terms. We then 
present a diagrammatic representation of the cumulant expansion that enables the derivation 
of the asymptotic convergence rate at any finite order of cumulant truncation. Much of the 
detailed examination of the convergence properties of the terms is given in two appendices. 
In Section III we review our key findings and discuss the expected numerical work required 
to implement the higher-order cumulants in actual calculations. 

II. THEORY 

A. Fourier path integral representations of the density matrix and partial aver- 
aging 

In this subsection we introduce the path integral representation of the quantum imaginary 
time density matrix, and demonstrate the utility of the Fourier representation of the paths 
in terms of a tail series. Some of the discussion in the subsection can be found in previous 
literature,'^*' but we find it useful to repeat the details for clarity and to establish the no- 
tation needed for the subsequent development. For simplicity we assume a one-dimensional 
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system with the extension to multidimensional systems to be presented with numerical ex- 
amples in a separate publication. 

We consider the Fourier representation of the Feynman-Kac formula for the quantum 
density matrixP 

p(x,x';P) = p fp (x,x';ft)K(x,x';ft) 

= Pf P {x, x'] ft) f[ ^= exp (~^ a fc) j ex P {-P V ( X ' x 'i f 3 )) (!) 

where pf p is the free-particle density matrix and ft = l/iksT) the inverse temperature with 
ks the Boltzmann constant. We have used the notation 

/= //[*(«)] du (2) 
Jo 

to represent an average of a function / with respect to the "imaginary time" variable u. In 
a similar manner, for a function g of two or more variables, we can write 

g= dui du 2 ... g[x{ux),y{u 2 ),..]. (3) 
Jo Jo 

Using this notation the time average of the potential in Eq.Q is 

l 

V(x, x', {a k }] ft) = duV x r (u) + era ■ A(u) , (4) 



with 



o 



x r (u) = x + (x' — x)u, (5) 



V m 

a= (ai,a 2 , ...), (7) 

A( M ) = (A 1 H,A 2 ( M ),...), (8) 

A fc (u) = V2 \ \ (9) 

and m is the particle's mass. The paths in Eq. Q, starting at coordinate x and ending at 
coordinate x' , are expanded in Fourier series, and each path is parametrized by an infinite 
set of Fourier coefficients (ai, a 2 , . . .). In order to calculate integral over the paths in Eq. ([T]), 
one needs to carry out an infinite dimensional integration over Fourier coefficients a±, a 2 , ■ ■ ■■ 
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Because an infinite dimensional integration is not numerically feasible, we truncate to a 
finite dimensional integral. The resulting finite dimensional integration converges to the 
exact result as the number of integration variables increases. To this end, we divide infinite 
dimensional integral Eq. into a finite dimensional integral over the first iV variables 
<2tv = (g&i, Q2, • • • j o/v) an d the tail integration (TI) over the infinite dimensional tail Sn+i = 

(Ojv+l, CJjV+2, • • •) 



K(x,x>;P) = f[(j ^g=exp (-± 




exp(-/?V)) 



(10) 



where the TI denotes the tail integration 



(exp(-/3V0) 



TI 



n 



k=N + l 




cxp 




x exp 



-(3 I <!u \ ' 

o 



x N (u) + aa N+1 ■ K N+ i(u) 



(11) 



In Eq.(ll) we have introduced the notation 

x N (u) = x r (u) + a N ■ A N (u) = x r (u) + a k^ki 



N 



(12) 



k=l 



The TI variables appear in the argument of potential function linearly, which can be ex- 
pressed as a scalar product 



a N+1 ■ A N+1 (u) = ^ a kA k {u). 



(13) 



k=N+l 



Equation (13) suggests a possible change of integration variables to a set of projection 
variables onto the vectors An+i(u), so that integration in the subspace orthogonal to this 
projection subspace might be performed trivially. The problem with a direct realization of 
this idea is that Eq. ([T]) contains an integration over u, and we have a continuum set of 
vectors Ajv+i(w). However, we can use the linearity of the integration operation allowing the 



order of integration in Eq.(ll) to be changed. To understand how the change of integration 



order helps evaluate Eq. (11), we expand the exponential function using a Taylor series 

<exp(-/3F)) T/ 



i + — 



p=l 



pi 



(14) 



where the TI of the p power of V defines the p order moment 

p. 




duV 



x N {u) + aa N+1 ■ A N+1 (u) 




dui - ■ ■ du p {T\V hc N {u k ) + cra N+1 ■ A N+ i(u k ) 

\k=l 



(15) 



TI 



In the second line of Eq. (15) we have interchanged the order of integration over u- variables 



and the TI variables ajv+i. At a fixed set of external time variables {ui, . . . , u p }, we now 
have a fixed set of vectors {A N+1 (ui), . . . , A N+ i(u p )} on which the vector a/v+i is projected. 
To proceed further, we recognize that in general, the A vectors are not orthonormalized. It 
is useful if we normalize the vectors and introduce non-orthogonal unit vectors 



9k 



= A N+1 (u k )/\/e(u k ), A; = 1,2, 



,P 



(16) 



where square of the normalization factor 



e(u k ) = 7(wfc,w fc ) = A N+ i(u k ) ■ A N+1 (u k ) = A^(w fc ) 

n=N+l 



(17) 



is a natural small parameter useful for much of the analysis found in the remainder of this 
work. We write 



e{u k ) 



1 

iV' 



(18) 



which, as explained and used in reference 1201 is a shorthand notation for asymptotic behavior 
of an integral with respect to u of the product of e(u) and a smooth function f(u) ; i.e. 

i 



due(u)f(u) = O 



N 



(19) 



The set of vectors {g k } defined in Eq. (16) is normalized but not orthogonal. It is useful 



to work with an orthogonal set that can by obtained from {g k } using Gramm-Schmidt 
orthogonalization procedure 

k k 



9k 



E 



a^ei, 22 a li = 1 > ( 20 ) 

i=i i=i 

where the vectors {e*i, • • • , e p } are expanded to ensure orthogonality; i.e. ej • e k = 5i k . The 
coefficients a ki and the vectors e*j can be easily calculated using the standard recurrence 



relations, so that we obtain 



fc-i \ 

^2a ki ei J /a k k, (21) 
i=i J 



&k = I <//,■ 

where 

olu = 9k -e h i = 1, • • • , k - 1, (22) 



\ 



k-l 



1 " • el) 2 (23) 

i=i 

Using these relations, we obtain, e.g., for the first three vectors 



where 



and 



e*i = gi, e 2 = (g 2 - g%\9\)l y 1 - g 2 i, 
e 3 = (<?3 - #3i<7i - (gs ■ 02)92)/ \j l ~9li- (93 • e 2 ) 2 , 
9s ■ e 2 = (#32 - 921931)/ yl -^21- ( 24 ) 



0i* =9i-9k= , llk (25) 



7*fc = l{ui,u k ) = A N+1 (ui) ■ A N+1 (u k ) 

N 

= min(«j, u fc ) - Mi?z fc - A n (ui)A n (u k ). (26) 

71=1 



Using the orthonormal set of constructed vectors, Eq. ( 13 ) becomes 

k 



a N +i ■ A N+1 (u k ) = ^/duk)^2a ki a N+ i ■ e { = y/ e(u k ) }^ a H ji (27) 



i=i t=i 
where denotes a projection of vector otv+i onto e«. 



The TI in Eq. (15) can be geometrically interpreted as an integration in an infinite- 
dimensional vector space whose elements are all possible vectors a^+i- If we choose 
{ei,-- - , e p } as the first p basis vectors, then an arbitrary vector will have the first p 
components {£1, • • • , Because the arguments of the potential function do not depend 
on the orthogonal projections, the integral with respect to the infinite set of components 
{£ p +i, £ p +2) " " " } m the directions orthogonal to the span of {e*i, • • • , e p } can be evaluated 
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analytically in Eq. (11) to give unity. As a result the TI in Eq. (15) is reduced to a 



p-dimensional Gaussian transform of the product of p potential functions 



G p (ui, ...,u p ) = J "J d u i - ■ ■ dUpG p (ui, . . . ,u p ) 



(28) 



o o 



G p (ui, . . . ,u p ) = ( Y[ y(x N (u k ) + aa N+ i ■ A N+ i(u k )) 

\k=l I 27 

P 



X 



n 

3=1 
V 

k=l 



X N [U k 



+ ay/e{u k ) 22akS 



i=i 



(29) 



We can modify the arguments of the potential functions in Eq. (29 ) to produce the partial 
average expression^ for the density matrix derived previously using alternate methods. We 
define A k using the expressions 



x N 



(u k ) + ay/ e(u k ) > y a H ii = x N (u k ) + ay/e(u k )£ k + A fe , 



(30) 



i=l 



where 

0, k = 1 

Yn=i a k£i + ("fcfc - l)6b, ^ = 2 
From the definition of the coefficients a k i, % = 1, . . . , k — 1, it is evident the coefficients are 



A fc = ay/e(u k ) 



(31) 



small for large N. Then for iV — > oo, we observe that at k ^ 2 



/fc-i fc-i 
A fc ~ ay/e(u k ) ^2 ~ ^2 



a 



ki 



(32) 



,i=i 



We can obtain the explicit asymptotic behavior for Eqs. (22 ) and (32 ) by using the following 
asymptotic expressions derived in Appendix A. 

1 



i i 



h{k) 




du 1 du 2 [l(u 1 ,u 2 )] k fi(u 1 )f 2 (u 2 ) 







O 
O 



k is odd 



k is even 



(33) 



where fi(ui), i = 1,2 are smooth functions and k = 1,2.... From Eqs. (22), (32) and 



asymptotic estimates Eq. (33), it follows that 



a k i ~ g k i = i(uk,Ui)/y/e{u k )e{ui), 

JL 2 JL 

^ ~ iV 3 ' ^ ~ TV 3 



and 



At 



7 2 K, Ml )/ V / ^FR ~ i/iv 3/2 . 



If we neglect the A& terms at k > 2 in (29) we find that the integration over £1, . . . ,£ p 
variables in Eq. (29) is reduced to a product of p one-dimensional integrals so that fi p = /xf. 



Consequently, the infinite power series is summed to produce an exponential function 

l + £^^M? = «p(-/?A*i), (34) 



P =i 



pi 



giving the so-called partial averaging (PA) formula first obtained by Doll et alP^ 



B. The cumulant expansion of the density matrix 



As shown elsewhere^ the partial average method expressed in Eq.(34) is the first term 



in the cumulant expansion of the density matrix. In this subsection we review the cumulant 
expansion and derive the asymptotic convergence characteristics of the various cumulant 
terms for the Fourier representation of the quantum density matrix. Much of what appears 
in the initial part of this subsection, in particular the defining relations for the cumulants, 
can be found elsewhere in the literature.^ However, we find it useful to spell out some 
well-known details to make the subsequent notation and discussion clear. 



The cumulant expansion can be obtained from Eq.(|14|) by writing 

<exp(-/3F)) T/ = 



z2 ~^~v P = ex p<x 



P =i 



p\ 



(35) 



where 




(36) 



k=l 



The last line of Eq. (36) is obtained by collecting the terms with the same power k of the 



potential function [or parameter /3] in 



Hck 



E 



r=l 



\r+l 



E 



PlV" ,Pr£C k , 



k\ 

Pi! ■■■p r \ 



(37) 
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The summation of the integer indices p±, ■ ■ ■ ,p r ^ 1 is restricted by the relation 



O., = { (pi, ■ ■ ■ ,p r ) : ^2,Vj = k 

3=1 



(38) 



Using Eqs. (37) and (38), the relations between the first four cumulants and the moments 
are 

2 

I^c2 — ^2 — A 4 ! , 

/i c3 = /i 3 - 3/^2 Aii + 2/xf . 

Ai c4 = A*4 - (4^3^! + 3/i|) + 12//2/X? - 6yuf 



(39) 



Given that the PA approximation is the first-order cumulant, the errors in the partial average 
method are determined by the second, third and higher-order cumulant terms. In a similar 
fashion the error in the second-order cumulant is determined by the third, fourth and higher 



order terms. Using Eqs. (39), one can also express the moments in terms of the cumulants 



A*3 



A^ci, 

A*c2 + Mel > 

A*c3 + 3At C 2A t ci + Mci- 



(40) 



/i 4 = Ai C 4 + 4/i c3 /i c i + 3a4 + 6a* C 2A«c1 + /4i 

These equations relating the cumulants and the moments can also be interpreted as 
recurrence relations that allow the expression of the higher-order cumulants in terms of 
lower-order ones. A general expression for the p th order moment can be derived by comparing 
coefficients in the power series expansion 



1 + —\»p 



P =i 



p\ 



OO ju 



OO h. 

X 



j=0 J \k=l 



(41) 
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From Eq. (41) we obtain 

dP 



ftp 



j=i j \k=\ 

1! ^ fci! dxP 



i=0 



fcl=l 



91 X 



2=0 fci,fc 2 =l 



( p x k 1 +k 2 



dxP 



x=0 



pl 



fei,i)2=l 

fcl+fc2=P 



+ 



1 F 



r! z — ' fcx! . . . fc r ! 
fci ,...,fe r =i 

ki + ...+k r =p 



f^cki ■ ■ ■ ^ck r + ...+//; 



cl 



The general term in Eq. (42) contains the factor 

p\ 



(42) 



fcl!...fc r F 

which defines the number of ways of grouping p objects into r groups of sizes fci, k^, . . . , k r , 
when the order within each group does not matter. These combinatorics imply that Eq. 



(42) can be represented as a sum over all partitions of {1, . . . ,p} vertices 



r=l S , i,S 2 ,...,5 r 



(43) 



where we have introduced a partition Si, S2, ■ ■ ■ , S r of a set of natural numbers {1, . . . , p} . 
As discussed elsewhere,^ if a partition Si, S 2 , ■ ■ ■ , S r contains, respectively, fci, fc 2 , . . . , fc r 
elements such that fci + k 2 + . . . + k r = p, then /i c [Si] = ^ c kn ■ ■ ■ i^c [S r ] — Hck r - The term 



in Eq. (43) for r = 1 corresponds to fi cp . The resulting expression is 



r=2 Si,S 2 ,...,Sr 



(44) 



where the summation is performed over the "proper partitions" when r ^ 2. To illustrate 
the idea of partitions, we show in Fig. [TJ all possible partitions for case p = 4 by connecting 
vertices 1, . . . ,4 by lines. Figure [I] can be compared to the last line of Eq. (40). The first 
line from the top in Fig. [I] displays the "improper partition" Si = (1, 2, 3,4) at r = 1, which 



corresponds to the term /i C 4 in Eq. (40). The partitions Si, 5*2 at r = 2 divide vertices into 



two groups. The first group shown in the second line consists of a cluster of three vertices 
and a cluster containing a single vertex. The second group shown in the third line consists 
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of two clusters each containing two vertices. The fourth and fifth lines show the partitions 
Si, S2, S3 at r = 3 consisting of a cluster of two vertices and two clusters each containing a 
single vertex. Finally, in the last line we find the partition Si, S 2 , S 3 , S 4 at r = 4 consisting 
of 4 clusters each containing a single vertex. 



C. Asymptotic convergence characteristics of the cumulant terms 

A key concern of any numerical path integral approach is the rate of convergence to the 
exact result with respect to the number of path variables included. In previous work^^S we 
have examined the Fourier path integral method along with partial averaging and gradient 
partial averaging, and examined the asymptotic convergence rates with respect to the num- 
ber of included Fourier coefficients. As mentioned previously, the partial average method 
is the first-order term in the cumulant expansion. The motivation of the current work is 
to derive the general convergence characteristics of the cumulants truncated at any finite 
order. 

We ignore the important and interesting question of the convergence of the cumulant 
expansion itself when a finite number of Fourier coefficients are included. We make the 
assumption that if we truncate the cumulant expansion at a finite order, the resulting ex- 
pression for the quantum density matrix will converge to the exact quantum density matrix 
in the limit of an infinite set of Fourier coefficients. 

In the following we first examine the convergence of the first three cumulants, and then 
derive an expression for the asymptotic convergence rates for cumulants of arbitrary, finite 
order. The analysis that follows is complex, mainly because the number of terms in high or- 
der cumulants is large. Many of the terms in each order have mixed asymptotic convergence 
characteristics along with internal cancellations within each order. We find a diagrammatic 
approach simplifies the complexity of the problem. 

Because 1) we are interested in the asymptotic convergence rate with respect to the 
number of included path variables N, 2) e(u) ~ 1/A, and 3) the arguments of the potential 



functions in Eq. (29) contain e(u), we make use of the Taylor series expansion for V around 
the point XN(iik)- For simplicity we assume that the potential function V has a convergent 
Taylor series to all orders. As a minimum requirement the potential energy must be infinitely 
differentiable. 
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1. The first-order cumulant 



Using the Taylor series expansion, we have 

1 oo 



du 



cxp 



V 



i 

„ oo 

J du Y. 



k=0 



v a 2k 



y {k KMu)) ak£k/2{u) 



£ 

A:=0 



(2Jfe)!! 



duV {2k) (x N (u))e k (u), 



(45) 



where the first line in Eq.(45) contains the "partial averaged potential," Vpa [x jv(u)], which 
is the Gaussian transform of the ordinary potential 



V pa {xn{u)} 



d£ 



exp 



e 



V 



Xn[U 



and where we have denoted V^ k \x) = d k V{x)/dx k and (2k)\\ = 2-4 



2k 



(46) 



2 k k\ if 



= 1,2,... [(2/c)!! = 1 at k = 0]. The Gaussian integral in Eq. (45) can be evaluated 
analytically resulting in the expression 



<<; exp 



2tt 



1 



if k is odd 
{k — 1)!! if k is even 



(47) 



Here (2Jfe-l)!! = 1 • 3 • . . . • (2k - 1) if k = 1, 2, . . . [(2Jfe-l)!! = 1 at k = 0]. Because integrals 
of e k {u) multiplied by a smooth function results in terms of order 1/N k , the last line in Eq. 



(45) gives the sought asymptotic expansion in powers of 1/N. 

There is an alternate derivation of the asymptotic expansion that in later developments 
we find to be particularly convenient for further generalizations. We rewrite the potential 
function as 

V x N (u) + o\Je(u)i 

where d = a— is a differential operator. Because <i-operators generate a commutative 

dx 

algebra similar to that of c-numbers, the resulting operator integral over £ can be calculated 



exp ^Je{u)^d) V(x N (u)) 



(48) 
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as an ordinary Gaussian integral, 



00 2 
9i{u,d) = J ^!= ex P + 
—00 

= exp f • (49) 



Applying the right-hand-side of Eq. (49) to the potential function, we obtain 



1 1 
Mci = / dugi(u,d)V(xN(u)) = f duexp ( ^—^-d 2 ) V(xn(u)). (50) 



Using a Taylor series to expand the exponential, the result is the asymptotic expansion for 
the PA potential 



V pa (xn(u)) = exp ( d 2 ) V(x N (u)) 



E " \, V(xn(u)), (51) 

fc=0 



formally identical to Eq. (45). The primitive Fourier path integral method corresponds to 
the first term, k = 0, while the "gradient partial average" approximation^^ is given by the 



first two terms (k = 0, 1) in the expansion of Eq.(51 ). 

Numerical evaluation of the one-dimensional w-integrals can employ any convenient 
quadrature rule. We have found the Gauss-Legendre quadrature formulaP' to be generi- 
cally useful 

N q 

A^ci = V PA (x N (u)) = ^WiVp A {x N {ui)) (52) 

i=i 

with Wi being a Gauss-Legendre weight and N q the number of quadrature points. The 
evaluation of the integral Eq. (52) requires N q calls to calculate the potential function Vpa, 
whereas the tOj coefficients do not depend on Vpa and can be precalculated. 



2. The second- order cumulant 



The second-order cumulant includes both the average potential and the average of the 
square of the potential. The average potential has already been discussed in Section [II C 1 , 
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so we consider the average of the square of the potential function V, Eq. (29) at p = 2. 



Using the operator representation (48) twice, we obtain 



G 2 {u 1 ,u 2 ) 



V 



x N (ui) + ay/e(ui)£i V x N {u 2 ) + cr\ / e(u 2 )(a 2 i^i + a 22 & 



ri 



g 2 (u 1 ,d 1 ,u 2 ,d 2 )V(l)V{2) 



(53) 



where 



g 2 (ui,di,u 2 ,d 2 ) 




d£id£ 2 



exp 



2tt "V 2 
x exp ( y/slZidi) exp ( </e^(a 21 £i + a 22 ^ 2 )d 2 



(54) 



In Eqs. (53) and (54) we have introduced the notation 



V(i) = V(x N (ui)), e, L = e(ui), di = a—, i = 1, 2 



(55) 



The differential operators d% and d 2 act, correspondingly, on potential functions V(l) and 
V{2). The operators d\ and d 2 commute and generate a commutative operator algebra. 



Evaluating the Gaussian integral Eq. (54), we obtain 

g 2 (ui, d l7 u 2 , d 2 ) = exp (- ^id\ + e 2 d\ + ^ 21 d\d 2 



exp 



9 'Yiididj 



(56) 



where 721 is defined in Eq. (26). In the second line, we have made use of the symmetry 



property 7^' = 7^ and the relation £j = e(ui) = j(ui,Ui) = 7^. Using the defining relation 



expressed in Eq. (49), Eq. (56) can also be rewritten as 



g 2 (ui,di, u 2 , d 2 ) = gi(ux, d^g^u^ d 2 )f 2 {T 2l ) 



(57) 



where 



hiXii) = exp(r 2 i) 

r 2 i = 72idi4- 



(58) 
(59) 
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From Eqs. (28), (50) and (53) we have 



/'2 

G 2 (ui,u 2 ) 



G 2 (ui,u 2 ), 

g 1 (u 1 ,d 1 )g 1 (u 2 ,d 2 )f 2 (T 21 )V(l)V(2) 
f2(r 21 )V PA (l)V PA (2), 



(60) 



i i 




du x du 2 g x {u x , d x )g x {u 2 , d 2 )V(l)V(2) 
J J du 1 du 2 V PA {l)V PA {2) 





1 1 



(61) 







and 



where 



Vc2 = Hz- n\ = G c2 (u l7 u 2 ) 
G c2 (u u u 2 ) = f c2 (T 21 )V PA (l)V PA (2) 

/ c2 (r 2 i) = / 2 (r 21 ) - 1 = exp (r 21 )-i. 



The integrand in \i\ is separable in the U\ and w 2 -time variables 

V PA {x N { Ul ))Vp A (x N (u 2 )) = V PA {1)V PA {2). 



(62) 



(63) 



(64) 



The integrations with respect to the variables fif over u± and u 2 are independent, and we 
call such terms that factor uncorrelated. In contrast \i c2 is proportional to 721 

G c2 (u u u 2 ) = [exp(T 21 )-l)V PA (l)Vp A (2) 



k=l 

00 „k 



£|f d\V PA {l) d\V PA {2) 



(65) 



and 



721 



2 E 

n=N+l 

00 



fc=l 



smlvmM! ) sinl7mM 2 



7m 



E 

n=7V+l 



cos7m(w 2 — Mi) — cos7rn(-u 2 + «i) 



7TO 



(66) 
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cannot be factored into a product of two separate functions of the time variables U\ and 
«2- The function 721 establishes a time correlation between time variables U\ and U2- The 



integrand Eq. (65) does not contain uncorrelated terms. All the uncorrelated terms cancel; 
i.e. if 721 — > 0, then G C 2{u\,U2) —> and fi C 2 — > 0. 

The demonstrated dependence of /i C 2 on the operator T 2 i can be graphically illustrated 
using a diagrammatic notation that proves to be especially valuable in simplifying the algebra 
for the higher-order cumulants. For /i c2 the diagrams are shown in Fig. [2] (a). We let small 
solid circles represent vertices 1 and 2 with the corresponding vertex functions Vp^l) and 
Vpa(2). The line that connects vertices 1 and 2 in the middle diagram of Fig. [2] (a) represents 
the interaction via T2i, which acts on the vertex functions via the differential operators d\ 
and d%. The two lines connecting vertices 1 and 2 shown in the right diagram, correspond 



to T^i entering in Eq. (65) with the weight coefficient 1/2!. In general, the k th power of a 
r 12 term, with weight coefficients l/k\, can be conveniently depicted by a diagram with the 
k lines connecting the two vertices. Summing the diagrams with the number of connecting 
lines from k = 1 up to 00 produces a diagrammatic representation of the series expansion 



given in Eq. (65). Corresponding to these two- vertex connected diagrams, we define a 
second-order cumulant function /i C 2(r2i), with the characteristic property that ^2^21) = 
if T2i = 0. In other words, the cumulant function takes a zero value on the corresponding 
disconnected vertex diagram, which is defined as two vertices not connected by a line that 
corresponds to the limiting case T 2 i — > [the left diagram in Fig. [2] (a)]. 



Using the expansion given in Eq. (65), the second cumulant function can be written as 
a sum 

00 

/i c2 = Y] f$ (67) 



k=l 



where jx^ is the corresponding contribution from the k th order derivative of the potential 



function Vpa 



1 1 



(fc) 

Vc2 




du\du2 



k\ 



d\V PA {\) d k 2 V PA (2) 



(68) 







Expressions for the integrals given in Eq. (68) are derived in Appendix A [see Eq. (Al) 
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Using the results of Appendix A, we find 

2 k a 2k 



(k) 
/ i c2 



E 



k\ ^ (irni) 2 • • • (nni,) 2 

m,- ,n k =N+l K u K K ' 



X 



-i 2 



du sin(7TOiM) • • • sin(7rn/ c 'U 



d k V PA (x N {u)) 
dx k 



(69) 



from which it follows that each term ^ ^ and, as a consequence, their sum, Eq. (67), is 
non-negative. 



As discussed for the first-order cumulant [see Eq. (52 )], fi^ can be calculated numerically 
using the Gauss-Legendre quadrature formula 



(fe) 

A*c2 



N q 

E 

i,j =1 



(k) d k Vp A {x N {ui)) d k V PA {x N (u j )) 



dx k 



dx k 



(70) 



where the symmetric semi-positive definite matrix 



zu 



(fc) 



a 



2k 



k\ 



(71) 



does not depend on the potential function, and its matrix elements can be precalculated. 
The evaluation of the double integral sum in Eq. (70) requires N q calls to calculate the k th 



order derivative of the potential function Vp A . Using the Cholesky decompositio: 



(72) 



where rjW is a unique lower triangular matrix with positive diagonal entries, Eq. (70) can 
be rewritten as 



(fe) 

Vc2 



E 

i=l 



From Eq. (|33j) follows that ^ an d of the same order of magnitude, 1/A^ 3 , and 



E»i 

.3=1 
,( 2 ) 



(fc) d k V PA (x N (uj)) 



n 2 



dx h 



(73) 



„ - ,,W _i_ ,,( 2 ) _i_ n 

Vc2 — H>c2 + ^c2 + U 



N 5 



(74) 



Using the asymptotic estimates, Eqs. (A11)-(A15) and (A16)-(A22), derived in Appendix 
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A, the asymptotically leading terms ^ an d be written as 



i 1 



u (1) 




a 
o o 

oo 

S ( 

n=7V+l 

3vr 4 iV 3 



dUidU2 r )21 



dV PA (x N (ui)) dV PA (x N (u 2 )) 



dx 



2a' 



Tin) 



dV(x) _ dV(x f ) 



dx 



dV(x) 
dx 



+ 



+ 



2(-l) J V 2 dV(x') 



dx dx 



dV{x') 
dx 



+ 



dx 

2 



dx 



+ 



1 

i\T5 



1 

AT5 



(75) 



and 



i i 



(2) 



a 

Y 




2 rf 2 \/ PA (a;Ar(Mi)) d 2 Vp A (x N (u 2 )) 
dU\dU2^2\~ 



dx' 



dx' 



o o 



a 



67T 4 iV 3 



(j 2 VpA(Xjv( w)) 

dx 2 



n 2 



1 

AT5 



(76) 



Combining the above asymptotic formulas, we obtain 



3vr 4 iV 3 



dx 



+ 



dx 



+ 



a 



du 



d 2 V PA (x N ( u)) 
dx 2 



n 2 



2(-l) J V 2 dK(x') 



1 

AT5 



(77) 



7T 4 7V 4 da; cte 

In the first line, we have collected terms of third order, while in the second line the terms of 



fourth and higher orders are included. Equation (|77l) confirms previous workP^ that shows 



that the PA approximation has a third-order convergence rate. Our numerical investigation 
of the PA convergence rate in a ID modeP^ is consistent with this asymptotic convergence 
rate. The calculation of the convergence constants in Eq. fl77| ) (the expression in the curly 
brackets) can be reduced to a one-dimensional integration over u. 



Equation ( 70 ) is comprised of N 2 terms composed of a product of k -order derivatives of 



the potential energy. The product of the two derivatives of the same order implies that for 
each k, the k th -oider derivative of Vp A [x^{ui)\ at each grid point Ui needs to be evaluated 



only once. In that way the N 2 operations required to evaluate ^ m Eq. (170b require 



only N q evaluations of the derivatives of the potential energy. Because the computational 
work can be expected to be dominated by the evaluation of the potential energy and its 
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derivatives, the work required to determine Eqs. ( 70 ) and ( 73 ) should scale nearly linearly 



in N q . From Eq. (|74|) only and required to attain an Af order asymptotic 



,(2) 



convergence rate further emphasizing that we can attain the same asymptotic convergence 
rate as the full cumulant with truncated algorithms that scale nearly linearly in N q . 

The above conclusions about the convergence rate are implicitly based on the assumption 
that the higher order cumulant terms, fi c k, k ^ 3, contribute at a faster rate. We check this 
assumption for the third cumulant term in the next section where we demonstrate that fi c3 
is fifth order in 1/N. 

3. The third-order cumulant 



Using methods identical to those employed in the derivation of Eq. (60), we obtain for 
the third-order moment 



A*3 = G z {u 1 ,u 2 ,u z ) 
i i i 




du 1 du 2 du 3 g 3 (ui, d x , u 2 , d 2 , u 3 , d 3 )V(l)V(2)V(3), 



ooo 



where 



oo oo oo 



93 




d |^#e X p(-^; + g +f |) + viT6i 



-oo — oo — oo 



+ V / ^( Q! 21^1 + "226)4 + V^^l^l + "326 + "33^)^ 



Evaluating the Gaussian integral Eq. (80), we obtain 



#3 = exp I - ^Teidf + Jijdidj j = exp ( - j 

\ i=i i<j=l / \ i,j=l / 

' 3 



where 



and 



/s(ri2, r23) r i3 ) 



/3(ri2,r 23 ,r 13 ) =exp ( ^ rv, 

\i<j=l 



(78) 
(79) 



(80) 



(81) 



(82) 



(83) 
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Using Eqs. (49) and (, 



81), we then obtain 



G 3 { Ul ,u 2 ,u 3 ) = f 3 (r 12 ,r 23 ,r 13 )v PA (i)v PA (2)v PA {3) 



(84) 



and 



where 



Vc3 = Hz - + 2f4 = G c3 (u u u 2 , u 3 ), 

G c3 ( Ul ,u 2 ,u 3 ) = fc3(Ti 2 ,T 23 ,T 13 )Vp A (l)Vp A (2)V PA (3) 



(85) 



U = / 3 (ri2, r 23 , r 13 ) - / 2 (r 12 ) - / 2 (r 23 ) - ; 2 (r 13 ) + 2. 



(86) 



We now show that the uncorrelated terms, — 3/i 2 /ii + 2[x\ sum to zero in the expression 
for /i c3 . We introduce the convenient notation 



Ti = r 23 , r 2 = r 13 , r 3 = r 12 . 



(87) 



and rewrite /, 



c3 



Ui^T^Ts) = exp X) r -5Z ex P(^) + 2 



,i=i 



i=l 



gj + r 2 + r s) fe - rj - r 2 fe - r* 



fc=i 



E 

k=2 



fc! 

( \ 

pPlpP2pP3 -pk 

1 1 1 2 1 3 l^i=l 1 i 



E 



Pl,P2,P3^0 
\P1+P2+P3=fc 



Pl!p 2 !P3! 



fc! 



(88) 



/ 



The inner summation in Eq. (88) is restricted to non-negative integer indices P\,p 2 ,p 3 with 



the additional requirement Pi + p 2 + Ps = fc. It is clear from the second line of Eq. (88) that 



the term at = 1 is zero. If any two of three indices take a zero value, say, p\ = p 2 = 
then p 3 = k and the corresponding term 

pPl pP2 pP3 



P X \p 2 \p 3 \ 



(pi=0,p2=0,p3=fc) 



1 3 

fc! 



(89) 



cancels in Eq. (88). We then find 



/- = E E 



fc=2 Pl,P2,P3&C k 



pPl pP2 pP3 

Pi!p 2 !p 3 ! 



(90) 
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where the following constraints on the indices are imposed 

. Pl,P2,P3>0, Pl+P2+P3 = k, 
C k = (P1,P2,P3) = \ (91) 

at least two of indices are non-zero. 



As a consequence of the constraints, we find Eq. (90) is proportional to the cross terms only 
and does not contain terms arising from the uncorrelated terms, — 3//2/-*i + , the terms 
with at least two zero indices. 

As with the second-order cumulant it is useful to introduce a diagrammatic notation for 
the various contributions. We define a diagram with three vertices 1, 2, 3 and corresponding 
vertex functions. In Fig. 2 (b), the three vertices are shown to be connected by lines 



that correspond to the pairwise vertex couplings T 12 , r 23 , T 31 . As follows from Eq. (86) or 



(90) the third cumulant function /x c3 (r 12 , r 2 3, r 31 ) is zero whenever the underlying diagram is 
disconnected. We find disconnected diagrams when any two of the vertex couplings are zero. 
For example, we obtain zero if Ti 2 = r 23 = 0, from which it follows that /i c3 (0, 0, T 3 i) = 0. 
In this case vertex 2 is disconnected from vertices 1 and 3, which are connected by a line 
that corresponds to the T 3 i ^ coupling. 

We now examine the cross terms up to third order; i.e. the terms with k = 2, 3 in Eq. 



(90). We find 



° pPl pP2 pP3 

k=2 P i, P2 ,p 3 ec k F1 ^ Fd 

= r^a + r 2 r 3 + + r 1 r 2 r 3 
+ l - (r 2 r 2 + v^l + r 2 r 3 + r 2 r 2 + r 2 r 3 + r^ 2 ) . (92) 

In Fig. [2] (b), we display the diagrams that correspond to the terms in the second line of 
Eq. ( |92| ). Using the asymptotic formulas derived in Appendix A [Eqs. ( |A31 ) and (A32) 



we find that terms in the second line of Eq. (92) are fifth order 

Y 1 T 2 + T 2 T z + T l T z + T 1 T 2 T z 

= 723713^1^2^3 + 713712^1^4 + 723712^1^3 

+7l27237l3^1<^3 5 ( 93 ) 

and the terms in the third line are of the higher order; namely, they produce a sixth order 



contribution. For example, if we examine a typical term of the third line in Eq. (92) 

r?r 2 = 7 2 2 37i3^24 (94) 
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and considering Eqs. (33), we find that 



7f 3 7i3 



1/N 6 



and 



r^r 2 ~ 



1 

7V6' 



(95) 



(96) 



After substituting Eq. (93) into (92) and (85), we find 



i i i 



V>c3 



a 




du 1 du 2 du 3 712713723 



d 2 V PA {\) d 2 V PA {2) d 2 V PA {3) 



dx 2 



dx 2 



dx' 



+3a 





1 1 1 

4 




duidu 2 du 3 7i 2 7i3 



d 2 V PA (l) dV PA (2) dV PA (3) 



dx 2 



dx 



dx 



+ 



1 

iV6 



(97) 







In Eq.(97), the first integral is a contribution from the loop diagram, the second line of Fig. 



2] (b), corresponding to the T^l^rsi term, whereas the second one represents a combined 
contribution, which is included by a factor 3 in front of the integral, from the degenerate 
diagrams shown in the first line of Fig. [2] (b) . 



As with the second-order cumulant terms, Eq. (97) implies that we can attain a N 
asymptotic convergence rate by retaining only the first two terms on the right-hand side. The 
work required when the integrals are evaluated numerically using Gauss-Legendre quadrature 
with N q quadrature points for a single one-dimensional imaginary-time integration, results 
in a total of 2N g calls of the first and second derivatives of the potential function Vp A . 



From the asymptotic estimates Eqs. (A31), (A32) obtained in Appendix A for the above 



fifth-order contributions, the integrals expressed in Eq. (97) can be rewritten as 

i 



yU C 3 




+ 3 



d 2 V(x) 
dx 2 



dV(x) 
dx 



2 



+ 3 



d 2 V(x') 
dx 2 



dV(x') 
dx 



(98) 



The calculation of the convergence constant in Eq. (98) is seen to be reduced to a one- 
dimensional integration over u. 



D. Higher-order cumulants and the linked-cluster theorem 



In principle, to calculate or estimate the asymptotic behavior of the higher-order cumu- 
lants fi c k, k > 4, one can use Eq. (37) as used in the previous sections to examine /i c2 and /i c3 . 
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However, with increasing cumulant order the number of terms to be estimated in Eq. (37) 



grows rapidly. Moreover, as we found in previous sections there are certain rules associated 
with connected diagrams, which help identify the terms that give non-zero contributions to 
the cumulants. The purpose of the present section is to generalize the results obtained for 
the second and third-order cumulants to higher orders. In particular, we derive asymptotic 
estimates for the higher order cumulants. 

By inspection of the results obtained for the first three moments fi p , p = 1,2,3, Eqs. 



(50), (56), (60), (78), and (81), we can generalize the integral expressions for the moments 



to an arbitrary index p: 



l l 



duf- du p g p {ui, d 1 , ■■■ ,u p , d p )V(l) ■ ■ ■ V{p) 



(99) 



o o 



where 

g p (ux,d!,--- ,u p ,dp) 



oo oo 



-oo — oo 



d£ d£ ( 1 P v k \ 

(2^/2 P ex P ( - 2 E $ + E e h*** S J 

exp (\ Yl nAd^j ■ ( 



Equation (100) is derived in Appendix B. 



It is useful to examine the diagonal terms in the exponent of Eq. ( 100 ) 

p 



i=l 



Y[ L(ui,£)v(i) =Y[v PA ( 



(101) 



4 = 1 



so that we can rewrite Eq. ( 99 ) as 



i i 



V P = j ' ■■■ J dui---du p f p {{Tij} p i<j= ^j V PA (l) ■ ■ ■ V PA {p) 



(102) 



o o 



where 



(103) 



\i<i=l 



The total number of couplings T k = T^, where i < j and i, j = 1, . . . ,p, p > 2, between all 
possible pairs of p vertices is 



K P = C p 



2 _ PiP ~ l ) 



(104) 
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We now define L\ and L 2 to be an arbitrary partition among p vertices such that the first 
cluster contains k± vertices and the second one k 2 vertices, with the total number of vertices 
being k\+k 2 = p. Without loss of generality, we assume that the first cluster contains vertices 
with numbers from 1 to k±, L\ = {1, . . . , ki}, and the second cluster contains vertices from 



k\ + 1 to p, L 2 = {ki + 1, . . . ,p}. Using this partition we call the exponent of Eq. (103) the 
"Hamiltonian" H\ 2 of the system of p vertices (using an analogy with a real Hamiltonian 
system), which can be divided as 

v 

Hi 2 = ^2 Fi i 

i<j=i 

= H 1 + H 2 + V 12 (105) 

where 

ki 

i<j=l 

V 

h 2 = T v ( 106 ) 

i<j=k±+l 

are, respectively, the "Hamiltonians" of the first and second clusters. The expression 

Vi2= ^ ( 10? ) 

is the "interaction" between the clusters. 

In accord with other diagrammatic approaches in physics, we define a diagram that has 
all vertices connected to be linked. If L\ and L 2 are two linked diagrams, we call L\ + L 2 
linked if the resulting diagram is linked. In contrast, if the interaction between the two 
diagrams V\ 2 = in Li + L 2 , the resulting diagram must be disconnected or unlinked. In 



the unlinked case, the factorization property given in Eq. (108) is valid 



tip = fikxVk 2 - (108) 



Additionally, if Eq.(108) is satisfied, then diagrams corresponding to the product of two 
moments must be unlinked. 

In previous sections we have shown that the contributions to the second and third-order 
cumulants are zero when the corresponding diagrams are disconnected. The cancellation of 
disconnected diagrams found in the second and third cumulants is an example of the linked 
cluster theorem, which applies not only in the present context, but in quantum field theory, 
many-body theory and cumulant expansions in statistical physics. 
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1. The linked- cluster theorem 



We now prove that only the terms that can be represented by linked diagrams make 
non-zero contributions to the corresponding cumulants. This statement is often called the 
linked-cluster theorem. We have already checked directly this statement for the second 
and third-order cumulants. Using the diagrammatic approach discussed in the previous 
subsection, we are now in a position to prove the linked cluster theorem for an arbitrary order 
cumulant term using the method of mathematical induction. We assume that cumulants 
from the second up to the p th order take zero values on the corresponding disconnected 
diagrams. Further, without loss of generality we assume that a disconnected p+1 cluster 
of vertices is made of two linked clusters L\ and L2 of sizes k% and &2, respectively, such 
that hi + hi — p + 1. Consequently, the interaction between the two clusters V12 = 0. It is 
straightforward to generalize the proof to the case of a disconnected cluster that is made of 
several linked clusters. In a second induction step, we must prove that the same property 
remains valid for this disconnected, p+l-vertex diagram; i.e. we must prove that ^ c (p+i) = 0. 



Rewriting Eq. (44) for the cumulant term of order (p + 1), we obtain 



p+1 

r=2 Si,S 2 ,...,S r 



(109) 



In view of the factorization property Eq. (108), we have 



(110) 



The summation in Eq. (109) 



p+1 

E E 

r=2 Si,S2,...,S T 

covers all possible "proper" partitions Si, S2, ■ ■ ■ , S r , r > 2 among p+1 vertices. Some 
members of these partitions, that we identify using the notation Sdi SC G {1, . . . ,p+l}, contain 
vertices from both clusters, L\ = {1, . . . , ki} and L 2 = {ki + 1, . . . ,p + 1}. Consequently, 
diagrams corresponding to these partition members are disconnected. Moreover, because 
r > 2 the maximum size of Sdisc, the maximum possible number of vertices in Sdi SC , is less or 
equal to p, and by induction, we have \i c [Sdisc] — 0. We then find that non-zero contributions 



to the sum in Eq. ( 109 ) are expected to arise only from partition clusters that belong either 



27 



to L\ or L2. Equation (109) can be rewritten as 



fea 



Mc(p+1) - Vkx^k 2 - E E A'c S 1 ! 



(1) 



il=l 5(1) oC 1 ) 
1 lj 

x E E ^ 

i a =l 5(2) ... 5(2) 
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(2) 



. . . fM c 



fJL c 



(1) 



5, 



(111) 



where S 1 ^, . . . , and , . . . , are, repectively, partitions of Li and L 2 . Finally, from 
the representation given in Eq.(43), it follows that A* c (p+i) = 0. The linked-cluster theorem 
is then proved. 



7 (i) 



7 (2) 



7(2) 



Asymptotic convergence rates for the truncated cumulant expansion 



According to the linked-cluster theorem, the expression given in Eq. (37) for cumulants 



includes only the terms that correspond to linked or connected diagrams. For this reason, 



we must omit terms in (37) having r > 2, the terms made of products of moments, whose 
diagrams are disconnected. The p th order cumulant can then be written as 

1 1 

A*cp = j ■■ J du!---du p f cp ({r fc }f4 j Vpa(1) • • • V PA {p) (112) 







where 



ftp = exp 



fc=0 




(r? 



Pi ysPK p 



■ ■ ■ V 



= E E — ^ ( 113 ) 

k=p-i Pl +-+p Kp =k F1 

Here, the subscrit c equivalently denotes "cumulant" and "connected". The inner sum- 
mation over non-negative integers Pi,P2, ■ ■ ■ ,Pk p is performed under the restrictions that 
Pi + P2 + • • • + Pk v — k and that all the vertices in the expansion should be connected by 
Tk = Tij couplings. The connectness of vertices is stressed by the subscript c. Using the 
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last restriction, the outer summation over k, starts at p — 1, because p vertices cannot be 
connected by less than p — 1 pairwise lines or couplings. 

The lowest order contribution to Eq. ( |113[ ) is the term k = p — 1 

U= E rw-iVi (ii4) 

linked diagrams 

where r^,^,..., _ a are pairwise couplings such that all p vertices are linked. Among all 
possible vertex connections we distinguish consecutive and centered ones. Figure [2] (c) illus- 
trates the notion of consecutive, centered, and loop diagrams in case of p = 4. Consecutive 
connections are those whose vertices, ordered in an arbitrary way, are linked consecutively 
one after another. If we label a set of ordered vertices from 1 to p, a consecutive connection 
corresponds to the following chain of operators 

r 12 r 23 • • • r(p_i) p . (115) 

If a vertex is linked to all other p — 1 vertices, then such a connection will be called centered. 
The corresponding "centered to vertex 1" operator chain takes the form 

r 12 r 13 ---r lp (116) 

For a consecutively linked cluster [Eq. ( |115 )], the contribution to the cumulant term is 



1 i 

cons 2(p— 1) 



tip™ = ° > ■■■ din-'-dup 712723 • • • 7(p-i) P 







x dV PA (l) d 2 V PA (2) d 2 V PA (3) d 2 V PA (p - 1) dV PA (p) 
doc doc dioc doc doc 



The numerical evaluation of the p-dimensional integral Eq.(117) with the Gauss-Legendre 
quadrature formula requires in total 2N q calls to calculate the first and second derivatives 
of the potential function Vp A , where N q is the number of quadrature points in a single 
dimension. 



In Appendix A, we derive an asymptotic estimate [Eq. (A44)] for integrals of the type 



Eq. (A33), from which we obtain an estimate for 



(Jcons 

,,cons P /iio\ 

»cp ~ jp^i (US) 

where C^ ons is a convergence constant. Moreover, in Appendix A we estimate the contribu- 
tion [Eq. ( |A45 )] from the loop diagram 

ri2r23 ■ ■ ■ ^ {p-i)jX pi-, (119) 
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which includes p T-couplings. We find that the corresponding loop diagram contributes 

1 1 

dui--- du p 712723 ■ ■ ■ 7(p-i) P 7pi 



= ° 2p 







n 



d 2 V PA {i) 
dx 2 



(120) 



to the cumulant. According to Eq.(A48), Eq. (120) is of the same order as the consecutive 
term; i.e. 



r^loop 
loop P 



CP J\T2p-l 



121) 



For a centered diagram representing Eq. (116) we obtain 



i i 



^cent = a 2(p-l) 



"cp 



J " J dui---du p 712713 • • • 7i P 







X 



d^VpAjl) 

dx'f- 1 



n 

i=2 



dVp A (i) 



dx 



(122) 



In Appendix A, we estimate this type of integral [(A51)] finding 



cent 



cent ^ P 



(123) 



Collecting all estimates Eqs. (118), (121), and (123) we can conclude that the total p th 
cumulant term scales as 



f^cp 



AT2 P -l 



(124) 



From Eq.(124), one can obtain an asymptotic estimate for the convergence rate of the 
cumulant expansion. If the cumulant expansion is truncated at order p, then the convergence 
rate of the truncated cumulant expansion is defined by the asymptotic behavior of the next 
cumulant term of order (p + 1). The resulting term scales as N~^ 2p+1 \ a central result of 
this work. 



III. DISCUSSION 

We have developed an asymptotic analysis of the cumulant expansion for the Fourier 
path integral representation of the quantum, imaginary-time density matrix. Starting from 
the Feynman-Kac formula expressed in Eq. ([l]), we have used the perturbation series given 
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in Eqs. (10) and (14), expanded in powers of the path averaged potential energy defined in 
Eq. (15). For the moments fi p , p = 1,2, .. . generated by this perturbation expansion we 



have found the compact integral representation expressed in Eq. (102). The integrand of the 



p-dimensional integral given in Eq. (102) is defined as the exponential function [Eq.(103) 
of a linear combination of differential operators acting on the product of p potential 
functions Vp^(l), • • • Vpa(p)- The operator T^- is proportional to the 7^ function [Eq.(26)], 
which asymptotically behaves as 7^ ~ 1/A^ 3 , which in the limit iV — > 00 goes to 0. In 



the same limit the exponential operator function given in Eq. (103) reduces to the identity 
operator. Using the relation e ~ l/N, we can conclude that Vpa, defined in Eq. (51), 
tends to V. Consequently, in the limit of large N we obtain a factorization of /x p — > [V] P , 
which results in Taylor's series expansion for the original exponential representation given 



in Eq. (|14j). 

When calculating the Feynman-Kac formula, it would be impractical to use the per- 
turbation series truncated at p — p max , because \i p considered as a function of iV does 
not go to zero as iV — > 00. As we have demonstrated, the cumulant expansion given in 



Eq. (36) behaves well asymptotically. In particular the p -order cumulant term fx cp scales 
as jV - ^ 2 ^ 1 ). At large enough N, the truncated cumulant expansion is expected to provide 
a good approximation to the exact density matrix, with the error scaling as j\M 2pmax+1 ). 



Using the polynomial expansion in products rf 1 ■ ■ ■ T P ^ P of powers of T k = T t j = jijclid 



for the exponential operator function given in Eq. (113), the integrand of the p-dimensional 



integral given in Eq. (112) can be represented as a sum of products of p potential functions 
Vpa(1), ■ ■ ■ Vpa(p) and their derivatives times the corresponding polynomial of 7-functions. 
Only the linked diagrams corresponding to the chains of operators I^ 1 • • • T P ^ P contribute 
to the final result. A numerical evaluation of the p-dimensional integral can be performed 
using, e.g. a p-dimensional Gauss-Legendre quadrature formula, with N q quadrature points. 

To find an optimal numerical scheme to evaluate the cumulants, it is important to estimate 
the amount of numerical work required. A reasonable estimate of this time is the number 



of potential function calls N ca u required to compute the integral given in Eq. (112). To 
minimize the time required to calculate the potential energy function at a fixed argument, 
it is useful to have an analytic expression for the Vpa function. It is well known that 



the Gaussian transform expressed in Eq. (46) can be evaluated analytically either for a 
polynomial or a Gaussian-type potential function. In general, we can assume that V can 
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be fitted by a finite combination of polynomial and/or Gaussian-type potential functions. 
For example, we have shown that a fit to the Lennard- Jones potential using two Gaussians 
gives numerical results that are indistinguishable to within statistical fluctuations of Monte 
Carlo path integral simulations^. 

As suggested in the previous paragraph, the computational cost in path integral simula- 
tions is dominated by the computational overhead required to evaluate the potential energy. 
Efficiency gains resulting from improved asymptotic convergence rates must be balanced 
with any possible increase in the number of calls required to evaluate the system potential. 



We have shown in Eq. (112) that we can attain the asymptotic convergence behavior at a 
given cumulant order p without including all terms in the series expansion for the cumu- 
lant. By truncating the expansion we are able to minimize the number of potential energy 
evaluations without sacrificing the improved asymptotic convergence rate. An important 
example is the case of the second-order cumulant discussed in the main text. We can attain 
the N~ 5 asymptotic convergence rate by including only the two terms represented in Eqs. 



(70) and (74). By virtue of the truncation, N ca u = 3N„, resulting in a scaling that is linear 



in N q . As a result of Eq. (112) this linear scaling in N q is maintained at all cumulant orders 
p so that N ca u = (r max + l)N q where 

'"max is the maximum order of the derivative of the 
potential function Vpa needed in the expansion to ensure the asymptotic convergence rate is 
jy-(2p+i) Even though the current paper has examined the simple case of one-dimensional 
systems, by using the expansion, the remarkable linear scaling in N q is maintained even for 
many-particle systems. Of course, if N q must be increased with cumulant order to obtain 
sufficiently accurate evaluations of the integrals with respect to the imaginary time variable, 
the scaling with cumulant order might be more severe than linear. Only numerical expe- 
rience will enable a full understanding of the scaling with N q . However, from the purely 
formal results, we expect to be able to extend partial averaging to higher order cumulants 
in a variety of important quantum problems. 

In reference [20j we have numerically investigated convergence characteristics for the en- 
ergy and heat capacity calculated with the first cumulant term in a one- dimensional Lennard- 
Jones model. Effects of the higher order cumulant terms and their convergence properties 
for iV-body systems are the subject of a separate publication. 
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Appendix A: Asymptotic estimates for integrals 

The purpose of this Appendix is to derive the asymptotic behavior of the integrals en- 
counted in the main text. First, we examine the behavior of the two-time integral defined 



in Eq.(33) in the asymptotic limit (N — > oo) 



1 1 

h(k) = J J du 1 du 2 {y(u 1 ,u 2 )] k fi{ui)f 2 {u 2 ) 



oo „i 



5^ 7 \5 7 \W dui sin(7T7iiWi) •••sin(7rn fc Mi)/i(ui) 



ni,— ,n k =N+l 

x / du 2 sin(7miM 2 ) ■ ■ ■ sm(nn k u 2 ) f 2 (u 2 ) (Al) 
Jo 

where fi(u) and f 2 {u) are smooth functions of u and k > 1 is an integer. In the second 



line we have used expansion given in Eq. (66) for 7. Expanding both fi and f 2 (denoted 



generally by f\ )2 ) in a Fourier series we obtain 

00 

/i,2(«i,2) = ^2 /i.2( m i,2) cos(7rmi i2 Mi l2 ). (A2) 



mi. 2=0 



If Eq.(A2) is substituted into Eq.(Al), the trigonometric integrals can be evaluated analyt- 



ically resulting in the expression 



n r --n fc =V+l v 17 mi,m 2 =0 

x J{m 1 ,n 1 , ■ ■ ■ ,n k )J(m 2 ,ni, • • • ,n k ) (A3) 
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where 



J{uq, ni, • • • , n&) = / du cos(irn Q u) si^nnxu) ■ ■ ■ sm(7rnku) 



1 



2 k +H k 



E 



O"! ■ • • CTfc 



1 if 1/ = 



0-0,0-1, — |Cfc=±l 



%TXV 



if z/ ^ 



(A4) 



i=o 



We assume that fi, i = 1, 2 are quadratically integrable. For such functions, to order e, 
these functions are well represented by a finite series with M £ \^ terms 

1/2 



r ( Mel - 2 ~ 

/ I /l,2(«) - E /l,2( m ) C0S 



irmu) 



^ e. 



(A5) 



In the asymptotic limit we recognize the condition iV > M e \^- Because the integral in Eq. 



(A4) is real, we observe the separate contributions for even and odd values of k. At k = 1 



and 2 we obtain, respectively, 



J(n ,rn) 



I _ ^_\ym+ni 
(vrni) (1 - (n /ni) 2 ) 



7rni 



E 

j=0 



7li 



and 



1 



J{nQf Til, ^2) {^ni,n.2+Wo ~^ ^m,7i2— no "I - ^n2,ni+no "I - $ri2,ni— no} 



(A6) 



(A7) 



where <5„ jm is the Kronecker delta. After substituting Eq. (A6) into (A3), we obtain 

-l)ji+h 



Mi) = £ 



■ 2 ■ 



n=7V+l 



7m 



E 



Jl J2=0 



_ 1) n / (2i 1 ) (l) _ / ^ (0) 



(2ji), 



[(-i)V 2 (2j2) (i)-/f 2) (o) 



(A8) 



where the following equalities have been used 

Me 



m=0 



£(-l) m /(m)m 2 ' 

E/( 



mm 



2j 



r-iv 



/<*)(!) 



H (2i) (o). 



m=0 



7T 



2.7 



(A9) 
(A10) 
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The functions on the right-hand side of Eqs. (A9) and (A10) are £-approximants using a 



truncated Fourier series. The leading term in Eq. (|A8|) is seen to be of third order 



E 

n=N+l 



(Tm) 



K-i)Vi(i)-/i(o)] 



x [(_i)" /2 (l)-/ 2 (0)]+O(iV- 5 ). 



(All) 



Equation (All) can be rewritten 
h(l) 



where 



1 



[/i(o)/ 2 (o) + A(i)/ 2 (i)] 



vr 4 [3N 3 

S 4 (N) [/!(l)/ 2 (0) + /i(0)/ 2 (l)]} + 0(N- 5 ), 



(A12) 



-11 



(A13) 

it,- 

n=N+l 

The asymptotic behavior in N can be estimated separately for even and odd values of N. 
At even N, we have 



oo .. oo 

S »W = E TlvT^I " E 



1 



1 



- f 

fc-i i L^i 



1/N 



E 



1/iV 



n=l ( 1 + 2 — J "=0 



1 + 



2n + 1 



N 



dx 



dx 

\J 1/N (l + 2x) k ' J (1 + i + 2x) fc 



1 



2iV fe! 



(A14) 



whereas at AT odd, one finds the same estimate by modulus, but with the opposite, "+" 
sign. Consequently, at any large A" we obtain 

_ 1)N +i 



S k (N) 



2N k 



(A15) 



After substituting Eq. (A7) into (|A3|), we obtain 

1 



J 2(2) = J2 



AL 



n=N+l 



(Tin) 



/i(0)/ 2 (0) + -^/ 1 (m)/ 2 



171} 



m=l 



oo 1 f _ _ 1 Me _ _ 

n=V+l ^ ' L m=l 



l-(^) 2 



(A16) 



771=1 
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where M £ = min(M £ i, M e2 ). After integrating by parts, we find 



K 4 



duf[ 2j \u)f 2 (u)= / duh{u)g j \u 



(A17) 



m=l 



where j — 0, 1, • • • . Using Eq.(A17), Eq. (A16) can be rewritten as 

I 2 (2) = K 7-^I + 3A 'i E r^6 + " 

(7m 4 (7rn) b 

n=iV+l v > n=N+l v y 



with 



and 



K = J duf 1 (u)f 2 (u), 
o 

i 

^=y ^a (2) h/ 2 (u) 





We observe that 



i i 




duidu 2 ^ 2 (ui,u 2 ) = 2j 



1 







n=iV+l 



7rn 



3tt 4 AT 3 ' 



(A18) 



(A19) 



(A20) 



(A21) 



so that the leading term of (A18) is of the third order, while the next term is of the fifth 



order. Moreover, as iV — > 00 the "normalized" function 7 2 (w 1 ,-u 2 ) function becomes a Dirac 
delta function 

l 2 (ui,u 2 ) 



lim — - 

J J duidu 2 ^ 2 {ui,u 2 ) 




5(ui - u 2 ). 



(A22) 



Equation (A22) has been proved in Appendix B of Reference fT9l 



For k > 2 the number of terms in Eq. (A4) grows exponentially with k, and the corre- 



sponding expressions for asymptotic constants quickly become quite unwieldy. However, we 
can still obtain asymptotic estimates at k ^ 3 by replacing the summations involved with 
integrations. By virtue of the Euler-MacLaurin summation formula, 35 these replacements 



become exact in the limit that iV — > 00. To demonstrate the approach, we examine Eq. (A3) 
noticing that the summation indices n\, • • • , begin with N + 1. We next replace these 
summation indices with continuous variables x\ = ri\/{N + 1), • • • ,Xk = Uk/(N + 1) and 
approximately replace the summation over indices by an integration over the correspond- 
ing x- variables with the integration range being defined from 1 to 00. For odd values of 
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k = 2p + 1, p = 0,1, ■ ■ ■ , the J-f unctions in Eq. (A3) are inversely proportional to A ~ N, 
so that asymptotically we obtain 

oo oo 

„ , , , If f dx-i ■ ■ ■ dxh , . Ch , . . 

I 2 (k = 2p + 1) ~ — J ■ ■ ■ J x \... xl k M*n -,x h ) = j^. (A23) 

11 

Here, there is no need to specify the function tp k ; the only requirement is the fmiteness 



of the integral. At even values of k — 2p + 2, p — 0, 1, ■ • • , the J-functions in Eq. (A3) 
are proportional to Kronecker deltas, which asymptotically tend to Dirac delta functions as 
N — > oo. As a result, we obtain k — 1 independent integration variables, which we set to be 
x 2 , ■ ■ ■ ,Xk, with x\ being a function of these k — 1 variables. The asymptotic expression for 
the integral takes the form 

oo oo 

_ . If f dxo ■ ■ ■ dxh , . Ch , , . 

I 2 (k = 2p + 2) ~ — J ■ ■ ■ J x l__ xl k M^ ■ ■ ■ ,x h ) = j^. (A24) 

11 

We next examine the asymptotic behavior of the three-time integrals of the type en- 



counted in Eq. (97) 



J 3 (l, 1,1) = J J J du 1 du 2 du 3 f(u 1 )f(u 2 )f(u 3 )>y(ui, u 2 )^{u 2 , u 3 )^{u x , u 3 ) (A25) 
o 

and 

i 

7 3 (0,1,1) = jjj du 1 du2du 3 f 1 (ui)f2(u 2 )f 2 (u3)j(ui, u 2 ) r y(u 1 , u 3 ) (A26) 
o 



Using the Fourier expansion of Eq. (A2), we obtain asymptotic results accurate to fifth 
order 

mi,m,2,m3=0 ni,n2,n3=7V+l 
X (&n2,ri3+m\ ~\~ ^n2,no y —m\){^n^ ) n\+m2 ^n^,n\— m^) (^«2,ni+m3 ^n2,n\—m^) 

' ' £ f(m 1 + m 2 )f(m 1 )f(m 2 ) + -\f(0)Y\ (A27) 



mi,m2=0 
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and 



73(0,1,1)= £ AW^Kl/^N E u^u^p 



»1 



7i {n\ — ml) 



n 2 



x [i _ _ o 

E 7 ^6 E A( m l)/2(^2)/2(m 3 ) 

*— ' (7rnr ' 

n=JV+l mi,rrt2,m3=0 



With some algebra it can be verified that 
i 



du [f(u)f 



mi,m2=0 



(A28) 



(A29) 



£ /lW/aH/alrns) [l - (-l)-^^] [l - (-1)"^] 

mi,m2,m3=0 

= A(o)[/ 2 (o)] 2 + /i(i)[/ 2 (i)] 2 - (-i) n / 2 (o)/ 2 (i)[A(o) + Mi)} (A30) 



Using Eqs. (A29) and (A30), Eqs. (A27) and (A28) can be rewritten 



/3(1,M) 



1 



du [f(u)f 



(A31) 



and 



/ 3 (o,i,i) 



[A(o)[/ 2 (o)] 2 + /!(i)[/ 2 (i)] 2 ] 



(A32) 



5ir 6 N 5 

We can also obtain an asymptotic estimate for the p-dimensional integral corresponding 
to a consecutively linked diagram 
1 1 

jcons = I ... I . . ■du p 'y{u 1 ,U2)'y{u2,U 3 ) ■ ■ •7(ttp_2,ttp_i)7(ttp-i,tt p ) 





X h{u 1 )f 2 {u 2 )h{Uz) ■ ■ ■ f 2 (u p - 2 )f2(Up-l)fl(u p ). 



(A33) 



Substituting the expansion given in Eq. (66) for the 7- function, we obtain 

Ji(ni) Ji(n p -i) 



rcons 
V 



= 2 t 

x J 2 (ni,n 2 )J 2 (n 2 ,n 3 )--- J 2 ( 



(A34) 
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where 



Ji(n) 



du sin(7rnu)/i(u), 



(A35) 



J 2 (ni,n 2 ) = / du [cos n (rii — n 2 )u — cosn (rii + n 2 )u]f 2 (u). 



(A36) 



Assuming the indices n,n\,n 2 > N + 1, A 7 " — > 00 are large, one can estimate the first 
integral using integration by parts. The second integral can be estimated by neglecting 
the contribution from the highly oscillating cos7r(ni + n 2 )u function. We then obtain the 
following estimates 

/x(0) - (-l)Vi(l) 



Ji(n) 



J 2 (n l7 n 2 ) ->■ J 2 (n 1 - n 2 ) = J du cos[7r(ni - n 2 )u}f 2 (u) 





(A37) 
(A38) 



Using the Fourier expansion [Eq. (A2)] and the property of orthogonality of the cosine 
functions, the last integral can be reduced to 



J 2 (n 1 - U 2 ) = ^[1 + £[ni-na[,o]/2(K - n 2 \) 



(A39) 



Substituting Eqs. ( A37[ ) and ( |A39[ ) into (jA34j) , we obtain 

2 3-p 
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V 



E 

rii,...,np_i=AT+l 



Lfi(o) - (-i) ni A(i)][/i(o) - (-1)^/1(1)] 



(7TOi) 3 (7m2) 2 



™ P -2) 2 (™ P -l) 3 



X [1 + 5[ m -n3|,o]/2(|rii - n 2 \) ■■■[1 + 5| np _ 2 _ np _3| i o]/ 2 (|np_2 - np-i|), (A40) 

It is convenient to replace the set of integer variables n 1; . . . ,n p _i running from N + 1 to 
infinity to the set of diagonal and off-diagonal integer variables 



(ni, n 2 , . . . , -)> (m, Ara 2 = n 2 - ni, . . . , An p _i = %-i - "1), 



so that the denominator in Eq. (A40) in the new variables becomes 



D = (7rni) 3 (7r[ni + An 2 ]) 2 ■ ■ ■ (tt[tii + Anp_ 2 ]) 2 (7r[ni + An p _ 1 ]) : 



When the diagonal variable n\ ^> |An 2 | 
and obtain to leading order 



An p _i|, we can set An 2 



D = (71-ni) 



2p 



(A41) 



(A42) 



Ara p _i = 



(A43) 
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The arguments of functions in the second line of Eq. (A40) depend only on the off-diagonal 



variables. For quadratically integrable functions the Fourier coefficients f 2 (\An\) — > as 
| An | — > oo, and to e error, the function can be well approximated by a finite Fourier series, 



Eq. (A5), with M £ + 1 terms. The second line of Eq. (A40) makes sizable contributions to 
the sum only at finite values of the off-diagonal variables |Ari2|, • • • , |An p _i|, restricted from 
above by the number M £ . This upper bound implies that asymptoticaly at 3> M £ the 



approximation for denominator, Eq. (A43), is justified and Eq. (|A40|) can estimated to give 



c E 



m=N+l 



■2p 



o 



1 



AT2p-l 



(A44) 



This asymptotic estimate is consistent with Eqs. ( A12 ) and ( A32 ), obtained in the particular 
cases p = 2 and 3. 

We can also examine the integral that corresponds to a loop diagram 



jloop 



dui ■ ■ ■ du p ^(u 1 ,u 2 ) r r(u2,u 3 ) ■ ■ -j(up-i, u p )^(u p ,ui) 



x /M •■■/(%) 



(A45) 



Substituting the expansion of the 7-function as in Eq. (A40), the loop integral can be 
reduced to 



jloop 

p 



E 



1 



(7ml) 2 • • • (irn r ^ 2 



m,...,n p =N+i v ±y v " ' v! 
x J(m - n p )J(n 1 - n 2 ) ■ ■ ■ J{n p ^ 2 - JV-0 J{ n p-i ~ n p) 



(A46) 



where 



./(;/) / du cos(irnu)f(u) = i [l + 5\ n \ fi ] f(n) 



(A47) 



Here, the f{n) 7 s are the Fourier coefficients in the Fourier expansion of the function f(u 
Using similar methods, we obtain the same asymptotic relations for 

1 



jloop _ Q 

P ^ V AT2P-1 



(A48) 



Similar arguments can be made to obtain the asymptotic behavior for consecutive and 



loop diagrams given in Eqs. (A44) and (A48). In the case of a consecutive diagram we 



have two end vertices 1 and p. Integration with respect to the variables associated with 
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the end vertices, u\ and u p , gives two extra factors in the denominator D of Eq. (A40), 
(nrii) 2 — > (nrii) 3 and (irn p _i) 2 — » (7rn p _i) 3 . These limits arise from the asymptotic estimate 



given in Eq. (A37). For vertices not at the ends, we obtain an extra 7(1*1, u p ) function 
coupling vertices 1 and p, which adds an extra factor (im p ) 2 to the denominator of Eq. 
(A46). Using the diagonal approximation n\ — . . 



n„-i or ni 



n p for consecutive 



or loop diagrams, we find that in both cases, the denominator has the same scaling as 
{nni) 2p , from which the same asymptotic expressions, Eqs. (A44) and (A48), follow. We 
have already faced the consecutive and loop diagrams for the particular cases p = 2 and 
3, when we have observed the same asymptotic behavior for the integrals ^(l) and ^2(2) [ 



compare Eqs. (A12) and (A18)], and for the integrals /3(C), 1, 1) and ^(l, 1, 1) [compare Eqs. 



(|A32f and ( |A3lj )]. 

In the case of a centered diagram we need to estimate the integral 

1 1 



rcent 



X/lOl) 



dui ■ ■ ■ du p u 2 )t(wi, u 3 ) ■ ■ ■ j(ui, Up) 
f[/a(«i) 



i=2 



(A49) 



Using the asymptotic expansion for the 7 function [Eq. (A37)] and the Fourier expansion 
for function f\, we obtain 

00 Vp— 1 
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/ 2 (o) - (-1)^/2(1) 
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^2 J(m,m, . . . ,7ip_i)/i(m) 



(A50) 



m=0 



where J(m, ni, . . . , n p _i) is defined in Eq. (A4). Equation (A50) is similar in structure to Eq. 
(A3), and we use the method of rescaled variables X\ — rii/(N+l), . . . , = n p _i/ (N+l) 
in analogy to the method used in deriving Eqs. (A23) and ( A24[ ). Separately examining the 
cases for even or odd values of p, we find the same asymptotic behaviors in both cases 



rcent 



o 
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i\T2p-l 



(A51) 
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Appendix B: Proof of Eq. fllOOh 



Taking the Gaussian integral Eq. (100) over £1, • • • , £ p , we obtain 



9 P = exp 



i=l 



-I 2^ 



k=i 



I 1 P P 



i=l k,k'=i 



We next write the exponent of Eq. (Bl) as 



ockiOtk'idkdk' 



i=l k,k'=i 
P 



9 Ckk'dkdk> 



k,k'=i 



where the coefficients are defined by 



min (k,k') 



_ 1/2 1/2 
Ckk' = £fr y a ki a k'i- 

i=l 



(Bl) 



(B2) 



(B3) 



We assume temporarily that k' > k. After the scalar multiplication of Eq. (21) by the 
vector gv, we obtain 

k 

" (B4) 



9k'k = 9k' ■ 9k — 2_j a ki a k'i 
i=l 



Evidently, if k > k' Eq. (B4) remains valid if k <->■ A;' so that 



min (k,k') 
9k'k = 9kk> = a ki a k'i 



(B5) 



i=l 



at any fc and k' . Consequently, we have 



1/2 1/2 _ 

Ckk' — £fc Pfefc' = 7fcfc' 



(B6) 



which proves Eq. (100). 
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Figure Captions 



1. Graphical representation of set partitions that provide the coefficients in the cumulant 
expansion. The top diagram represents the improper partition S = (1,2,3,4). The 
second line represents all ways of partitioning three elements into one group with the 
fourth element remaining. The third line shows all possible pair partitions, the fourth 
and fifth lines represent partitions containing a single pair, and the final line shows 
the last possibility consisting of four clusters each with a single vertex. Each of the 



lines correspond in order to the 5 terms in Eq.(40) with the corresponding product of 
cumulants and associated coefficients indicated. 

2. Diagrammatic representation of the cumulants. Each solid circle represents a vertex 
function Vp^(z) and each line represents an interaction (see text for definitions). In 
(a) the contributions to the second-order cumulant is shown, with the first diagram in 
(a) being unlinked and contributing 0. Contributions to the third-order cumulant are 
shown in (b), with the first three diagrams representing terms in the second line of Eq. 



(92) and the last diagram being an example of a loop diagram. In (c) three fourth- 
order diagrams are shown demonstrating consecutive, centered and loop diagrams 
respectively. 
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FIG. 2: Kunikeev et al. 
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